%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Marginal Effect for model 5: with unobs and no school dummies %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%% To study the marginal effect, consider an individual in a specific
%%%% group. So first, find a group
D=100;

bnew1=bnew;
gg=71;  %% Find the groups whose mean characteristics is closest to the population.

unobs=dat(gg).unobs;

X_f = dat(gg).xt;

Mg = dat2(gg).Mg;


n=size(X_f,1);
mr=n;

k=17;

bnew=bnew1;
W_f= dat(gg).wt;
wt=W_f;
    indp_b1f=zeros(mr,1); 
    indp_b2f=zeros(mr,1);
    deri_f=zeros(mr,1); 

    for j=1:D
    
        bbf=X_f*bnew(1:k+0)+wt*X_f(:,1:k)*bnew(k+0+1:2*k+0)+Mg(:,j)*bnew(2*k+0+1)+bnew(2*k+0+2)*ones(mr,1)+exp(bnew(2*k+0+3))*unobs(:,j)*ones(mr,1); %unobs is 1 by D for each group.         
        indp_b1f=indp_b1f+ 1./(1+exp(-2*bbf));
       
    end;
    temp3=indp_b1f/D;
    
    F=temp3;  % the F average over all 100 simulations.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  
 naivemarg=[2*diag(F.*(1-F))*(bnew(1)+.2*X_f(:,1)*bnew(17)) 2*F.*(1-F)*[bnew(2:16)' bnew(18+0:35+0)']];  
% not included age_2, constant as separate X here. But still include it in WX, but not report the results of its marginal effect (the 2nd to the last one).
naivemargm=mean(naivemarg(1:n,:))'*100;  

ini_n=zeros(n,k-1);
X_fa=X_f;  %%alias
for member=1:n  
   
    %%%%% Marginal effect for the discrete variables

    for kk=3:k-1   %%%to k-1 ,as the last regressor is squared-age
        
           Mg1=Mg;   
           
         X_fa(member,kk)=0;   
         indp_b1f=zeros(mr,1); 
         indp_b2f=zeros(mr,1);
         deri_f=zeros(mr,1); 
         
         for j=1:D
                
        bbf=X_fa*bnew(1:k+0)+wt*X_fa(:,1:k)*bnew(k+0+1:2*k+0)+Mg1(:,j)*bnew(2*k+0+1)+bnew(2*k+0+2)*ones(mr,1)+exp(bnew(2*k+0+3))*unobs(:,j)*ones(mr,1); %unobs is 1 by D for each group.
       
        indp_b1f=indp_b1f+ 1./(1+exp(-2*bbf));
        
        end;
        
        
    temp3=indp_b1f/D;
    
    F_before=temp3;  % the F average over all 100 simulations.

         X_fa(member,kk)=1;        
         indp_b1f=zeros(mr,1); 
         indp_b2f=zeros(mr,1);
         deri_f=zeros(mr,1); 
   
         for j=1:D
                
        bbf=X_fa*bnew(1:k+0)+wt*X_fa(:,1:k)*bnew(k+0+1:2*k+0)+Mg1(:,j)*bnew(2*k+0+1)+bnew(2*k+0+2)*ones(mr,1)+exp(bnew(2*k+0+3))*unobs(:,j)*ones(mr,1); %unobs is 1 by D for each group.
                       
        indp_b1f=indp_b1f+ 1./(1+exp(-2*bbf));
       
        end;
        
        
    temp3=indp_b1f/D;
    
    F_after=temp3;  % the F average over all 100 simulations.

    
        amarg=F_after-F_before;
      
        ini_n(member,kk)=amarg(member);
        
        X_fa=X_f;
    end
end

ini_n2=ini_n(:,3:k-1);
ma1_n=mean(ini_n2)'*100; %By multiply by 100, express in terms of percentage point
naivemargm2=naivemargm;
naivemargm2(3:k-1)=ma1_n; %%naive estimation of marginal effects


    indp_b1f=zeros(mr,1); 
    indp_b2f=zeros(mr,1);
    deri_f=zeros(mr,1); 

    for j=1:D
    
        bbf=X_f*bnew(1:k+0)+wt*X_f(:,1:k)*bnew(k+0+1:2*k+0)+Mg(:,j)*bnew(2*k+0+1)+bnew(2*k+0+2)*ones(mr,1)+exp(bnew(2*k+0+3))*unobs(:,j)*ones(mr,1); %unobs is 1 by D for each group.
                       
        indp_b1f=indp_b1f+ 1./(1+exp(-2*bbf));
        indp_b2f=indp_b2f+ 1./(1+exp(2*bbf));
        
        %gradient
        deri=4./(exp(2*bbf)+exp(-2*bbf)+2);
        
        deri_f=deri_f+deri;
        
        
    end

fderiv_tanh=deri_f/D;   % the fderiv_tanh average over all 100 simulations.

ini=zeros(n,k-1);
affect=zeros(n,k-1);
X_fa=X_f;  
for member=1:n   
    %%%%% Marginal effect for the continuous variables
    kk=2;  %%% years in school
    dXdx=zeros(size(X_f));
    dXdx(member,kk)=1;
   
    %DXdx for contextual variables (not include dummied in contextual).
    dXdxW=zeros(n,k);
    dXdxW(member,kk)=1;

    
    partial=inv(eye(n)-diag(fderiv_tanh)*bnew(2*k+0+1)*W_f)*diag(fderiv_tanh)*(dXdx*bnew(1:k+0)+W_f*dXdxW*bnew(k+0+1:2*k+0)); %%% Calculate dM/db for group g
    amarg=2*diag(F.*(1-F))*(dXdx*bnew(1:k+0)+W_f*dXdxW*bnew(k+0+1:2*k+0)+bnew(2*k+0+1)*W_f*partial); %% marginal effect for ALL members
    ini(member,kk)=amarg(member);
    amarg(member)=0;
    affect(:,kk)=affect(:,kk)+amarg;

    kk=1; %%% age
    dXdx=zeros(size(X_f));
    dXdx(member,kk)=1;
    dXdx(member,kk+16)=2/10*X_f(member,1);
    
    dXdxW=zeros(n,k);
    dXdxW(member,kk)=1;
    dXdxW(member,kk+16)=2/10*X_f(member,1);

    partial=inv(eye(n)-diag(fderiv_tanh)*bnew(2*k+0+1)*W_f)*diag(fderiv_tanh)*(dXdx*bnew(1:k+0)+W_f*dXdxW*bnew(k+0+1:2*k+0)); %%% Calculate dM/db for group g
    amarg=2*diag(F.*(1-F))*(dXdx*bnew(1:k+0)+W_f*dXdxW*bnew(k+0+1:2*k+0)+bnew(2*k+0+1)*W_f*partial); %% marginal effect for ALL members
    
    ini(member,kk)=amarg(member);
    amarg(member)=0;
    affect(:,kk)=affect(:,kk)+amarg;

    %%%%% Marginal effect for the discrete variables

    for kk=3:k-1   %%%to k-1 ,as the last regressor is squared-age

        X_fa(member,kk)=0;
  
           Mg1=zeros(mr,D);
        
            for j=1:D;          
                
                x0=ones(mr,1)/2;
                euclid=1;
                while euclid>1e-7
                   xx_new=tanh(X_fa*bnew(1:k+0)+W_f*X_fa(:,1:k)*bnew(k+0+1:2*k+0)+W_f*x0*bnew(2*k+0+1)+bnew(2*k+0+2)*ones(mr,1)+exp(bnew(2*k+0+3))*unobs(:,j)*ones(mr,1));
                   euclid=max(sum((xx_new-x0).*(xx_new-x0)));
                   x0=xx_new;
                end
                
            Mg1(:, j)=wt*xx_new;
            
            end
         indp_b1f=zeros(mr,1); 
         indp_b2f=zeros(mr,1);
         deri_f=zeros(mr,1); 
         
         for j=1:D
                
        bbf=X_fa*bnew(1:k+0)+wt*X_fa(:,1:k)*bnew(k+0+1:2*k+0)+Mg1(:,j)*bnew(2*k+0+1)+bnew(2*k+0+2)*ones(mr,1)+exp(bnew(2*k+0+3))*unobs(:,j)*ones(mr,1); %unobs is 1 by D for each group.
        indp_b1f=indp_b1f+ 1./(1+exp(-2*bbf));
 
        end;
        
        
    temp3=indp_b1f/D;
    
    F_before=temp3;  % the F average over all 100 simulations.

         X_fa(member,kk)=1;
          
           Mg1=zeros(mr,D);
      
            for j=1:D;          
              
                x0=ones(mr,1)/2;
                euclid=1;
                while euclid>1e-7
                   xx_new=tanh(X_fa*bnew(1:k+0)+W_f*X_fa(:,1:k)*bnew(k+0+1:2*k+0)+W_f*x0*bnew(2*k+0+1)+bnew(2*k+0+2)*ones(mr,1)+exp(bnew(2*k+0+3))*unobs(:,j)*ones(mr,1));
                   euclid=max(sum((xx_new-x0).*(xx_new-x0)));
                   x0=xx_new;
                end
                
            Mg1(:, j)=wt*xx_new;
            
            end
       
         indp_b1f=zeros(mr,1); 
         indp_b2f=zeros(mr,1);
         deri_f=zeros(mr,1); 
   
         for j=1:D
                
        bbf=X_fa*bnew(1:k+0)+wt*X_fa(:,1:k)*bnew(k+0+1:2*k+0)+Mg1(:,j)*bnew(2*k+0+1)+bnew(2*k+0+2)*ones(mr,1)+exp(bnew(2*k+0+3))*unobs(:,j)*ones(mr,1); %unobs is 1 by D for each group.
   
        indp_b1f=indp_b1f+ 1./(1+exp(-2*bbf));
       
        end;
      
    temp3=indp_b1f/D;
    
    F_after=temp3;  % the F average over all 100 simulations.

        amarg=F_after-F_before;
    
        ini(member,kk)=amarg(member);
        amarg(member)=0;
        affect(:,kk)=affect(:,kk)+amarg;

        X_fa=X_f;
    end
end

ma1=mean(ini);  %%% average marginal effect of initiators
ma2=mean(affect)/n;%%% average marginal effect of followers
ma=[ma1; ma2]'*100%%% By multiply by 100, express in terms of percentage point

[naivemargm naivemargm2];
naive=naivemargm2





